Membrane Thickness Analysis Tutorial¶
Table of Contents¶
- Introduction
- Installation and Setup
- Pipeline Overview
- Input Data and Configuration
- Running the Pipeline
- Command Line Interface
- Results Analysis
- Results Visualization
- Extract Coordinates of Matched Measurement Points
- Convert Thickness Data into Motive Lists
Introduction¶
This tutorial demonstrates how to run the memthick pipeline for measuring membrane thickness from cryo-electron tomography (cryo-ET) membrane segmentations, as well as examples of how to analyze and visualize the output data.
Pipeline Stages¶
Input Preparation:
- Perform instance membrane segmentations and curate labels. We recommend using MemBrain-seg [1] with connected components analysis.
- For separating or merging labels post-segmentation, we recommend using the 3D lasso tool in the MemBrain v2 napari plugin or the napari-segselect plugin [2].
- Prepare labels (as a config file or dictionary), segmentation, and tomogram files.
- The input data used in this tutorial is an instance segmentation of a Chlamydomonas reinhardtii tomogram [3], containing nuclear envelope (NE), inner and outer mitochondrial membranes (IMM and OMM), generated using MemBrain-seg [1]. Input data can be found here along with an example configuration file.
- The pipeline can also be run on semantic segmentations (i.e., all membrane instances labeled with ones).
Surface Extraction: Extract surface point coordinates and orientations from segmentations using marching cubes, refine orientations using weighted neighbor averaging, and separate bilayer surfaces using principal component analysis of normal vector orientations.
Thickness Measurement: Match points between surfaces with geometric constraints (maximum thickness and cone angle search). Measure thickness as Euclidean distance between matched points.
Intensity Profiling: Extract intensity profiles from the input tomogram along vectors connecting matched points and validate measurements using profile features (user-configurable, for details see Stage 3: Intensity Profile Analysis).
Quality Filtering: Remove uncertain measurements based on intensity profile features such as dual minima detection, central maximum presence, and signal-to-noise ratios.
Analysis & Visualization: Comprehensive statistical analysis and visualization of results with interactive Plotly plots.
References:¶
[1] L. Lamm, et al. MemBrain v2: an end-to-end tool for the analysis of membranes in cryo-electron tomography. bioRxiv. 2025. https://www.biorxiv.org/content/10.1101/2024.01.05.574336v2
[2] napari-segselect: Napari Plugin to clean membrain-seg segmentations
[3] R. Kelley, et al. Towards community-driven visual proteomics with large-scale cryo-electron tomography of Chlamydomonas reinhardtii. bioRxiv. 2024 https://www.biorxiv.org/content/10.1101/2024.12.28.630444v1
Installation and Setup¶
Dependencies¶
# Core dependencies (should be installed in cryoCAT environment)
pip install pyyaml tqdm plotly pandas
conda install numba -c conda-forge
# Optional: GPU acceleration (recommended for large datasets)
conda install cudatoolkit
In [1]:
%load_ext autoreload
%autoreload 2
import memthick_final as memthick # Main pipeline
import memthick_analyze_plot as maple # Analysis and plotting module
from IPython.display import IFrame
Pipeline Overview¶
Core Parameters Reference¶
Input Files¶
| Parameter | Description | Default | Notes |
|---|---|---|---|
segmentation_path |
Input segmentation MRC file | (required) | Must have correct voxel size metadata |
tomo_path |
Tomogram MRC file for intensity analysis | None | Required for intensity profiling |
config_path |
YAML file with membrane labels | None | Alternative to membrane_labels |
membrane_labels |
Dict mapping names to label values | {"membrane": 1} |
e.g., {"NE": 1, "ER": 2} |
output_dir |
Output directory | Same as input | Where results are saved |
Surface Processing¶
| Parameter | Description | Default | Recommended | Impact |
|---|---|---|---|---|
mesh_sampling |
Marching cubes step size | 1 | 1-2 | Higher = faster, less detail |
interpolate |
Add interpolated surface points | True | ✓ | Improves surface coverage |
interpolation_points |
Points between vertices | 1 | 1-2 | Higher = denser sampling |
refine_normals |
Improve normal vectors | True | ✓ | Essential for bilayer separation |
radius_hit |
Normal refinement radius (voxels) | 10.0 | 5-15 | 10 should work well |
flip_normals |
Orient normals inward | True | ✓ | For thickness measurements |
Thickness Measurement¶
| Parameter | Description | Default | Recommended | Impact |
|---|---|---|---|---|
max_thickness |
Maximum search distance (nm) | 8.0 | 6-10 | Membrane-dependent |
max_angle |
Cone search angle (degrees) | 1.0 | 1-3 | Lower = more precise |
direction |
Measurement direction | "1to2" | Try both | Should not affect results |
use_gpu |
GPU acceleration | True | ✓ if available | Impacts only the thickness measurement step |
Intensity Profiling¶
| Parameter | Description | Default | Recommended | Impact |
|---|---|---|---|---|
extract_intensity_profiles |
Enable intensity analysis | True | ✓ | Required for automated filtering of thickness results |
intensity_extension_voxels |
Profile extension distance | 10 | 8-15 | Context around measurement points |
intensity_extension_range |
Analysis window | (-10, 10) | Match extension_voxels | Range for feature detection (setting this too wide can be problematic in case of e.g. membrane-bound ribosomes) |
intensity_min_snr |
Minimum signal-to-noise ratio | 1.0 | 0.5-2.0 | Peak prominence above baseline noise |
intensity_margin_factor |
Position tolerance | 0.0 | 0.1-0.3 | How much further out of the paired points can the minima be detected (0.3 ≈ c.a. 1-1.5 voxels) |
intensity_require_both_minima |
Both minima required | True | ✓ | Must find two minima within the measurement range |
intensity_central_max_required |
Require central maximum | True | ✓ | Two minima separated by a central maximum |
Input Data and Configuration¶
Example Data Structure¶
project/
├── inputs/
│ ├── segmentation.mrc # Instance or semantic segmentation
│ ├── tomogram.mrc # Original tomogram (for intensity)
│ └── config.yml # Membrane label configuration
└── outputs/ # Generated results
Configuration File Format¶
# config.yml
NE: 1
OMM: 2
IMM: 3
ER: 4
Setup Paths¶
In [2]:
# Input paths
segmentation_path = 'inputs/2738_z200to340_segmented.mrc'
tomo_path = 'inputs/2738_z200to340.mrc' # For intensity profiling
config_path = 'inputs/config_2738.yml'
# Alternative: Direct label specification
membrane_labels = {
'NE': 1,
'OMM': 2,
'IMM': 3
}
# Output directory
output_dir = 'outputs/'
In [3]:
# GPU-accelerated full pipeline with intensity profiling
results = memthick.run_full_pipeline(
segmentation_path=segmentation_path,
output_dir=output_dir,
config_path=config_path, # or membrane_labels=membrane_labels
# Surface processing
interpolate=True,
interpolation_points=1,
refine_normals=True,
flip_normals=True,
# Thickness measurement
max_thickness=8.0,
max_angle=1.0,
direction='1to2',
use_gpu=True,
# Intensity profiling (recommended)
extract_intensity_profiles=True,
tomo_path=tomo_path,
intensity_extension_voxels=8,
intensity_extension_range=(-8, 8),
intensity_min_snr=0.1,
intensity_margin_factor=0.3, # 30% tolerance for minima positions extending outward from the paired point positions
intensity_central_max_required=True,
intensity_require_both_minima=True,
# Optional outputs
save_vertices_mrc=False,
save_thickness_mrc=False
)
2025-08-20 14:56:42,428 - Starting full membrane thickness analysis pipeline for inputs/2738_z200to340_segmented.mrc 2025-08-20 14:56:42,429 - Intensity profiling enabled with tomogram: inputs/2738_z200to340.mrc 2025-08-20 14:56:42,430 - Step 1: Processing membrane segmentation 2025-08-20 14:56:42,432 - Reading segmentation from inputs/2738_z200to340_segmented.mrc... 2025-08-20 14:56:42,435 - Voxel size: 0.7840 nm 2025-08-20 14:56:42,436 - Origin: (np.float32(0.0), np.float32(0.0), np.float32(0.0)) 2025-08-20 14:56:42,437 - Shape (ZYX): (141, 1024, 1024) 2025-08-20 14:56:42,438 - Processing NE (label 1) 2025-08-20 14:56:42,505 - Extracting surface points with step size 1... 2025-08-20 14:56:42,506 - Extracting surface points with step size 1... 2025-08-20 14:56:47,626 - Extracted 371023 surface points 2025-08-20 14:56:48,422 - Extracted 371023 surface points 2025-08-20 14:56:48,423 - Interpolating surface points (before: 371023 vertices) 2025-08-20 14:57:04,281 - Interpolation: 371023 → 498069 vertices 2025-08-20 14:57:04,823 - After point interpolation: 498069 vertices 2025-08-20 14:57:04,825 - Refining normals with weighted-average method... 2025-08-20 14:57:04,827 - Building KD-tree for neighbor searches... 2025-08-20 14:57:05,052 - Separating surfaces using marching cubes orientation... 2025-08-20 14:57:05,052 - Separating bilayer surfaces... 2025-08-20 14:57:05,065 - Surface 1: 246748 points 2025-08-20 14:57:05,078 - Surface 2: 251321 points 2025-08-20 14:57:05,090 - Orienting normals... 2025-08-20 14:57:05,091 - Computing normal orientations... 2025-08-20 14:57:05,789 - Building adjacency graph in batches of 2000... 100%|██████████| 250/250 [01:07<00:00, 3.71it/s] 2025-08-20 14:58:13,284 - Computing minimum spanning tree... 2025-08-20 14:58:13,287 - Propagating orientations... 2025-08-20 14:58:13,288 - Refining normals using weighted average of neighbor normals... 2025-08-20 14:58:13,289 - Refining normal directions using single-scale approach... 100%|██████████| 250/250 [00:40<00:00, 6.19it/s] 2025-08-20 14:58:53,649 - Separating bilayer surfaces... 2025-08-20 14:58:53,659 - Surface 1: 244604 points 2025-08-20 14:58:53,666 - Surface 2: 253465 points 2025-08-20 14:58:53,667 - Successfully separated surfaces 2025-08-20 14:58:53,668 - Successfully separated bilayer surfaces: 2025-08-20 14:58:53,669 - Surface 1: 244604 points 2025-08-20 14:58:53,669 - Surface 2: 253465 points 2025-08-20 14:58:53,670 - Updating vertex volume with interpolated points... Updating vertex volume: 100%|██████████| 498069/498069 [00:00<00:00, 546419.02it/s] 2025-08-20 14:58:54,622 - Final vertex count: 498069 2025-08-20 14:58:54,622 - Saving outputs for NE: 2025-08-20 14:58:54,623 - Vertices: 498069, Surface 1: 244604, Surface 2: 253465 2025-08-20 14:58:59,927 - Saved CSV with 498069 vertices to outputs/2738_z200to340_segmented_NE_vertices_normals.csv 2025-08-20 14:58:59,928 - Processing OMM (label 2) 2025-08-20 14:58:59,952 - Extracting surface points with step size 1... 2025-08-20 14:58:59,952 - Extracting surface points with step size 1... 2025-08-20 14:59:02,637 - Extracted 131984 surface points 2025-08-20 14:59:02,937 - Extracted 131984 surface points 2025-08-20 14:59:02,940 - Interpolating surface points (before: 131984 vertices) /home/ms/deglushk/projects/tm_membranes/bilayer_thickness/memthick_final.py:1288: RuntimeWarning: invalid value encountered in divide interp_normal /= np.linalg.norm(interp_normal) 2025-08-20 14:59:08,335 - Interpolation: 131984 → 174727 vertices 2025-08-20 14:59:08,508 - After point interpolation: 174727 vertices 2025-08-20 14:59:08,509 - Refining normals with weighted-average method... 2025-08-20 14:59:08,510 - Building KD-tree for neighbor searches... 2025-08-20 14:59:08,578 - Separating surfaces using marching cubes orientation... 2025-08-20 14:59:08,579 - Separating bilayer surfaces... 2025-08-20 14:59:08,583 - Surface 1: 84659 points 2025-08-20 14:59:08,586 - Surface 2: 90068 points 2025-08-20 14:59:08,586 - Orienting normals... 2025-08-20 14:59:08,587 - Computing normal orientations... 2025-08-20 14:59:08,863 - Building adjacency graph in batches of 2000... 100%|██████████| 88/88 [00:18<00:00, 4.65it/s] 2025-08-20 14:59:27,787 - Computing minimum spanning tree... 2025-08-20 14:59:27,789 - Propagating orientations... 2025-08-20 14:59:27,790 - Refining normals using weighted average of neighbor normals... 2025-08-20 14:59:27,790 - Refining normal directions using single-scale approach... 100%|██████████| 88/88 [00:13<00:00, 6.29it/s] 2025-08-20 14:59:41,787 - Separating bilayer surfaces... 2025-08-20 14:59:41,797 - Surface 1: 76865 points 2025-08-20 14:59:41,802 - Surface 2: 97862 points 2025-08-20 14:59:41,802 - Successfully separated surfaces 2025-08-20 14:59:41,803 - Successfully separated bilayer surfaces: 2025-08-20 14:59:41,803 - Surface 1: 76865 points 2025-08-20 14:59:41,804 - Surface 2: 97862 points 2025-08-20 14:59:41,804 - Updating vertex volume with interpolated points... Updating vertex volume: 100%|██████████| 174727/174727 [00:00<00:00, 521552.66it/s] 2025-08-20 14:59:42,169 - Final vertex count: 174727 2025-08-20 14:59:42,170 - Saving outputs for OMM: 2025-08-20 14:59:42,170 - Vertices: 174727, Surface 1: 76865, Surface 2: 97862 2025-08-20 14:59:44,029 - Saved CSV with 174727 vertices to outputs/2738_z200to340_segmented_OMM_vertices_normals.csv 2025-08-20 14:59:44,029 - Processing IMM (label 3) 2025-08-20 14:59:44,052 - Extracting surface points with step size 1... 2025-08-20 14:59:44,053 - Extracting surface points with step size 1... 2025-08-20 14:59:47,730 - Extracted 235660 surface points 2025-08-20 14:59:48,227 - Extracted 235660 surface points 2025-08-20 14:59:48,228 - Interpolating surface points (before: 235660 vertices) /home/ms/deglushk/projects/tm_membranes/bilayer_thickness/memthick_final.py:1288: RuntimeWarning: invalid value encountered in divide interp_normal /= np.linalg.norm(interp_normal) 2025-08-20 14:59:57,585 - Interpolation: 235660 → 314937 vertices 2025-08-20 14:59:57,922 - After point interpolation: 314937 vertices 2025-08-20 14:59:57,923 - Refining normals with weighted-average method... 2025-08-20 14:59:57,924 - Building KD-tree for neighbor searches... 2025-08-20 14:59:58,059 - Separating surfaces using marching cubes orientation... 2025-08-20 14:59:58,059 - Separating bilayer surfaces... 2025-08-20 14:59:58,070 - Surface 1: 155861 points 2025-08-20 14:59:58,082 - Surface 2: 159076 points 2025-08-20 14:59:58,089 - Orienting normals... 2025-08-20 14:59:58,093 - Computing normal orientations... 2025-08-20 14:59:58,569 - Building adjacency graph in batches of 2000... 100%|██████████| 158/158 [00:38<00:00, 4.11it/s] 2025-08-20 15:00:37,056 - Computing minimum spanning tree... 2025-08-20 15:00:37,058 - Propagating orientations... 2025-08-20 15:00:37,059 - Refining normals using weighted average of neighbor normals... 2025-08-20 15:00:37,060 - Refining normal directions using single-scale approach... 100%|██████████| 158/158 [00:25<00:00, 6.15it/s] 2025-08-20 15:01:02,756 - Separating bilayer surfaces... 2025-08-20 15:01:02,778 - Surface 1: 139731 points 2025-08-20 15:01:02,782 - Surface 2: 175206 points 2025-08-20 15:01:02,782 - Successfully separated surfaces 2025-08-20 15:01:02,783 - Successfully separated bilayer surfaces: 2025-08-20 15:01:02,784 - Surface 1: 139731 points 2025-08-20 15:01:02,785 - Surface 2: 175206 points 2025-08-20 15:01:02,785 - Updating vertex volume with interpolated points... Updating vertex volume: 100%|██████████| 314937/314937 [00:00<00:00, 553723.64it/s] 2025-08-20 15:01:03,386 - Final vertex count: 314937 2025-08-20 15:01:03,386 - Saving outputs for IMM: 2025-08-20 15:01:03,386 - Vertices: 314937, Surface 1: 139731, Surface 2: 175206 2025-08-20 15:01:06,742 - Saved CSV with 314937 vertices to outputs/2738_z200to340_segmented_IMM_vertices_normals.csv 2025-08-20 15:01:06,749 - Step 2: Measuring membrane thickness 2025-08-20 15:01:06,749 - Processing thickness for NE 2025-08-20 15:01:06,750 - Reading segmentation from inputs/2738_z200to340_segmented.mrc... 2025-08-20 15:01:06,751 - Voxel size: 0.7840 nm 2025-08-20 15:01:06,751 - Origin: (np.float32(0.0), np.float32(0.0), np.float32(0.0)) 2025-08-20 15:01:06,751 - Shape (ZYX): (141, 1024, 1024) 2025-08-20 15:01:06,752 - Loading data from outputs/2738_z200to340_segmented_NE_vertices_normals.csv... 2025-08-20 15:01:07,158 - CUDA GPU is available and will be used 2025-08-20 15:01:07,159 - Starting thickness measurement (1to2)... 2025-08-20 15:01:07,159 - Measuring thickness from surface 1 to surface 2... 2025-08-20 15:01:07,160 - Starting GPU thickness measurement with 498069 points... 2025-08-20 15:01:07,160 - Source points: 244604, Target points: 253465 2025-08-20 15:01:07,161 - Max thickness: 8.0 nm (10.20 voxels) 2025-08-20 15:01:07,161 - Max angle: 1.0 degrees 2025-08-20 15:01:07,304 - Launching CUDA kernel with 1946 blocks and 256 threads per block... 2025-08-20 15:01:07,566 - Retrieving results from GPU... 2025-08-20 15:01:08,346 - Processing matches on CPU to ensure one-to-one matching... 2025-08-20 15:01:14,199 - Found 221676 valid thickness measurements 2025-08-20 15:01:14,202 - Mean thickness: 5.66 nm 2025-08-20 15:01:14,205 - Min: 0.78 nm, Max: 8.00 nm 2025-08-20 15:01:14,206 - Processing completed in 7.05 seconds 2025-08-20 15:01:14,223 - Statistics saved to outputs/2738_z200to340_segmented_NE_thickness_stats.log 2025-08-20 15:01:14,224 - Creating optimized thickness CSV with only valid measurements... 2025-08-20 15:01:14,307 - Saving 221676 valid thickness measurements to outputs/2738_z200to340_segmented_NE_thickness.csv 2025-08-20 15:01:19,184 - Thickness statistics: 2025-08-20 15:01:19,185 - Mean: 5.658 nm 2025-08-20 15:01:19,187 - Std: 0.688 nm 2025-08-20 15:01:19,187 - Min: 0.784 nm 2025-08-20 15:01:19,188 - Max: 7.995 nm 2025-08-20 15:01:19,188 - Membrane thickness analysis complete! 2025-08-20 15:01:19,192 - Thickness analysis for NE completed successfully 2025-08-20 15:01:19,192 - Step 3: Processing intensity profiles for NE 2025-08-20 15:01:19,192 - Validating segmentation and tomogram compatibility... 2025-08-20 15:01:19,193 - Reading segmentation from inputs/2738_z200to340_segmented.mrc... 2025-08-20 15:01:19,194 - Voxel size: 0.7840 nm 2025-08-20 15:01:19,194 - Origin: (np.float32(0.0), np.float32(0.0), np.float32(0.0)) 2025-08-20 15:01:19,194 - Shape (ZYX): (141, 1024, 1024) 2025-08-20 15:01:19,195 - Reading tomogram from inputs/2738_z200to340.mrc... 2025-08-20 15:01:19,201 - Tomogram voxel size: 0.7840 nm 2025-08-20 15:01:19,201 - Tomogram shape (ZYX): (141, 1024, 1024) 2025-08-20 15:01:19,202 - ✓ Segmentation and tomogram are compatible 2025-08-20 15:01:19,202 - Dimensions: (141, 1024, 1024) 2025-08-20 15:01:19,202 - Voxel size: 0.7840 nm 2025-08-20 15:01:19,203 - ============================================================ 2025-08-20 15:01:19,203 - Processing: 2738_z200to340_segmented_NE_thickness.csv 2025-08-20 15:01:19,203 - Tomogram: 2738_z200to340.mrc 2025-08-20 15:01:19,204 - ============================================================ 2025-08-20 15:01:19,204 - Loading thickness data... 2025-08-20 15:01:19,531 - Loading tomogram... 2025-08-20 15:01:19,532 - Extracting intensity profiles...
Reading tomogram from inputs/2738_z200to340.mrc... Tomogram voxel size: 0.7840 nm Tomogram shape (ZYX): (141, 1024, 1024) Normalizing tomogram using method: zscore Extracting intensity profiles with extension of 8 voxels Including physical coordinates with voxel size: 0.7840 nm Processing 221676 valid coordinate pairs Extracting intensity values for 221676 profiles...
2025-08-20 15:01:28,543 - Extracted 221676 intensity profiles 2025-08-20 15:01:28,544 - Filtering intensity profiles... 2025-08-20 15:01:28,544 - === Statistics of the extracted intensity profiles ===
Successfully extracted 221676 intensity profiles
2025-08-20 15:01:59,872 - Step 1: Found 164795 candidate profiles 2025-08-20 15:01:59,873 - Median separation: 4.87 ± 0.88 voxels 2025-08-20 15:01:59,873 - Median prominence SNR: 4.1x 2025-08-20 15:01:59,874 - === Apply filtering criteria: (1) two minima (with user-defined SNR compared to baseline) separated by max; (2) the min are positioned between the two matched points from which thickness was measured === 2025-08-20 15:02:43,513 - Step 2: 200594/221676 profiles passed the filtering criteria (90.5%) 2025-08-20 15:02:43,889 - Saving results to: outputs/
=== Summary of intensity profiles post filtering === Total profiles analyzed: 221,676 Candidate profiles found after step 1: 164,795 Profiles that passed criteria: 200,594 (90.5%) Profiles excluded: 21,082 (9.5%) === Characteristics of intensity profiles === Median separation of minima: 4.87 ± 0.88 voxels Median SNR of minima: 4.1x (minima are 4.1x more pronounced than the baseline) === Applied filtering criteria (user-defined) === Required SNR of minima: 0.1x (minima must be 0.1x more pronounced than the baseline) Position filtering: Exact (no margin) Central max required: True Extension range from midpoint: (-8, 8) voxels === Reasons for intensity profile exclusion === minima_outside_extended_region: 18,106 (8.2%) insufficient_snr_min1: 2,923 (1.3%) no_central_maximum: 39 (0.0%) insufficient_snr_min0: 14 (0.0%) === Quality metrics of included profiles === Mean separation of minima: 4.91 ± 0.90 voxels Median separation of minima: 4.87, Range: 1.85-13.44 Mean SNR of minima: 4.1 ± 2.0x Median SNR of minima: 3.7x, Q25-Q75: 2.7-5.0x (Minima are on average 4.1x more pronounced than the baseline) Saved cleaned thickness DataFrame: outputs/2738_z200to340_segmented_NE_thickness_cleaned.csv Merging Pass 1 characterization feature data with 221676 profiles Available feature data entries: 221676 Profiles with feature data: 221676/221676 Profiles with minima positions: 164795/221676 Profiles with separation distance: 164795/221676 Profiles with SNR data: 164795/221676 Profiles with central max position: 0/221676 Saved original intensity profiles with features: outputs/2738_z200to340_segmented_NE_int_profiles.pkl Mapping 200594 cleaned profiles to 200594 passing indices Merging Pass 2 filtering feature data with 200594 profiles Available feature data entries: 221676 Profiles with feature data: 200594/200594 Profiles with minima positions: 200594/200594 Profiles with separation distance: 200594/200594 Profiles with SNR data: 200594/200594 Profiles with central max position: 200594/200594 Saved cleaned intensity profiles with features: outputs/2738_z200to340_segmented_NE_int_profiles_cleaned.pkl Saved statistics: outputs/2738_z200to340_segmented_NE_filtering_stats.txt
2025-08-20 15:03:14,119 - Intensity profiling for NE completed successfully 2025-08-20 15:03:14,121 - Profiles extracted: 221676 2025-08-20 15:03:14,122 - Profiles after filtering: 200594 2025-08-20 15:03:14,123 - Filter pass rate: 90.5% 2025-08-20 15:03:14,124 - Processing thickness for OMM 2025-08-20 15:03:14,125 - Reading segmentation from inputs/2738_z200to340_segmented.mrc... 2025-08-20 15:03:14,127 - Voxel size: 0.7840 nm 2025-08-20 15:03:14,128 - Origin: (np.float32(0.0), np.float32(0.0), np.float32(0.0)) 2025-08-20 15:03:14,129 - Shape (ZYX): (141, 1024, 1024) 2025-08-20 15:03:14,129 - Loading data from outputs/2738_z200to340_segmented_OMM_vertices_normals.csv... 2025-08-20 15:03:14,291 - CUDA GPU is available and will be used 2025-08-20 15:03:14,291 - Starting thickness measurement (1to2)... 2025-08-20 15:03:14,292 - Measuring thickness from surface 1 to surface 2... 2025-08-20 15:03:14,293 - Starting GPU thickness measurement with 174727 points... 2025-08-20 15:03:14,294 - Source points: 76865, Target points: 97862 2025-08-20 15:03:14,294 - Max thickness: 8.0 nm (10.20 voxels) 2025-08-20 15:03:14,295 - Max angle: 1.0 degrees 2025-08-20 15:03:14,300 - Launching CUDA kernel with 683 blocks and 256 threads per block... 2025-08-20 15:03:14,413 - Retrieving results from GPU... 2025-08-20 15:03:14,416 - Processing matches on CPU to ensure one-to-one matching... 2025-08-20 15:03:16,188 - Found 70793 valid thickness measurements 2025-08-20 15:03:16,189 - Mean thickness: 4.63 nm 2025-08-20 15:03:16,191 - Min: 0.78 nm, Max: 8.00 nm 2025-08-20 15:03:16,192 - Processing completed in 1.90 seconds 2025-08-20 15:03:16,199 - Statistics saved to outputs/2738_z200to340_segmented_OMM_thickness_stats.log 2025-08-20 15:03:16,200 - Creating optimized thickness CSV with only valid measurements... 2025-08-20 15:03:16,230 - Saving 70793 valid thickness measurements to outputs/2738_z200to340_segmented_OMM_thickness.csv 2025-08-20 15:03:17,802 - Thickness statistics: 2025-08-20 15:03:17,804 - Mean: 4.632 nm 2025-08-20 15:03:17,805 - Std: 0.768 nm 2025-08-20 15:03:17,806 - Min: 0.784 nm 2025-08-20 15:03:17,807 - Max: 7.995 nm 2025-08-20 15:03:17,807 - Membrane thickness analysis complete! 2025-08-20 15:03:17,808 - Thickness analysis for OMM completed successfully 2025-08-20 15:03:17,809 - Step 3: Processing intensity profiles for OMM 2025-08-20 15:03:17,809 - Validating segmentation and tomogram compatibility... 2025-08-20 15:03:17,810 - Reading segmentation from inputs/2738_z200to340_segmented.mrc... 2025-08-20 15:03:17,812 - Voxel size: 0.7840 nm 2025-08-20 15:03:17,812 - Origin: (np.float32(0.0), np.float32(0.0), np.float32(0.0)) 2025-08-20 15:03:17,813 - Shape (ZYX): (141, 1024, 1024) 2025-08-20 15:03:17,814 - Reading tomogram from inputs/2738_z200to340.mrc... 2025-08-20 15:03:17,823 - Tomogram voxel size: 0.7840 nm 2025-08-20 15:03:17,823 - Tomogram shape (ZYX): (141, 1024, 1024) 2025-08-20 15:03:17,824 - ✓ Segmentation and tomogram are compatible 2025-08-20 15:03:17,825 - Dimensions: (141, 1024, 1024) 2025-08-20 15:03:17,826 - Voxel size: 0.7840 nm 2025-08-20 15:03:17,872 - ============================================================ 2025-08-20 15:03:17,873 - Processing: 2738_z200to340_segmented_OMM_thickness.csv 2025-08-20 15:03:17,873 - Tomogram: 2738_z200to340.mrc 2025-08-20 15:03:17,874 - ============================================================ 2025-08-20 15:03:17,875 - Loading thickness data... 2025-08-20 15:03:17,981 - Loading tomogram... 2025-08-20 15:03:17,983 - Extracting intensity profiles...
Reading tomogram from inputs/2738_z200to340.mrc... Tomogram voxel size: 0.7840 nm Tomogram shape (ZYX): (141, 1024, 1024) Normalizing tomogram using method: zscore Extracting intensity profiles with extension of 8 voxels Including physical coordinates with voxel size: 0.7840 nm Processing 70793 valid coordinate pairs Extracting intensity values for 70793 profiles...
2025-08-20 15:03:21,297 - Extracted 70793 intensity profiles 2025-08-20 15:03:21,299 - Filtering intensity profiles... 2025-08-20 15:03:21,299 - === Statistics of the extracted intensity profiles ===
Successfully extracted 70793 intensity profiles
2025-08-20 15:03:30,353 - Step 1: Found 25887 candidate profiles 2025-08-20 15:03:30,354 - Median separation: 5.09 ± 1.45 voxels 2025-08-20 15:03:30,355 - Median prominence SNR: 3.8x 2025-08-20 15:03:30,356 - === Apply filtering criteria: (1) two minima (with user-defined SNR compared to baseline) separated by max; (2) the min are positioned between the two matched points from which thickness was measured === 2025-08-20 15:03:43,648 - Step 2: 36971/70793 profiles passed the filtering criteria (52.2%) 2025-08-20 15:03:43,740 - Saving results to: outputs/
=== Summary of intensity profiles post filtering === Total profiles analyzed: 70,793 Candidate profiles found after step 1: 25,887 Profiles that passed criteria: 36,971 (52.2%) Profiles excluded: 33,822 (47.8%) === Characteristics of intensity profiles === Median separation of minima: 5.09 ± 1.45 voxels Median SNR of minima: 3.8x (minima are 3.8x more pronounced than the baseline) === Applied filtering criteria (user-defined) === Required SNR of minima: 0.1x (minima must be 0.1x more pronounced than the baseline) Position filtering: Exact (no margin) Central max required: True Extension range from midpoint: (-8, 8) voxels === Reasons for intensity profile exclusion === minima_outside_extended_region: 27,250 (38.5%) insufficient_snr_min1: 6,352 (9.0%) no_central_maximum: 168 (0.2%) insufficient_snr_min0: 47 (0.1%) zero_baseline_noise: 5 (0.0%) === Quality metrics of included profiles === Mean separation of minima: 4.79 ± 1.26 voxels Median separation of minima: 4.67, Range: 1.85-14.77 Mean SNR of minima: 3.3 ± 2.1x Median SNR of minima: 2.9x, Q25-Q75: 2.0-4.1x (Minima are on average 3.3x more pronounced than the baseline) Saved cleaned thickness DataFrame: outputs/2738_z200to340_segmented_OMM_thickness_cleaned.csv Merging Pass 1 characterization feature data with 70793 profiles Available feature data entries: 70793 Profiles with feature data: 70793/70793 Profiles with minima positions: 25887/70793 Profiles with separation distance: 25887/70793 Profiles with SNR data: 25887/70793 Profiles with central max position: 0/70793 Saved original intensity profiles with features: outputs/2738_z200to340_segmented_OMM_int_profiles.pkl Mapping 36971 cleaned profiles to 36971 passing indices Merging Pass 2 filtering feature data with 36971 profiles Available feature data entries: 70793 Profiles with feature data: 36971/36971 Profiles with minima positions: 36971/36971 Profiles with separation distance: 36971/36971 Profiles with SNR data: 36971/36971 Profiles with central max position: 36971/36971
2025-08-20 15:03:50,322 - Intensity profiling for OMM completed successfully 2025-08-20 15:03:50,324 - Profiles extracted: 70793 2025-08-20 15:03:50,324 - Profiles after filtering: 36971 2025-08-20 15:03:50,325 - Filter pass rate: 52.2% 2025-08-20 15:03:50,326 - Processing thickness for IMM 2025-08-20 15:03:50,327 - Reading segmentation from inputs/2738_z200to340_segmented.mrc... 2025-08-20 15:03:50,330 - Voxel size: 0.7840 nm 2025-08-20 15:03:50,330 - Origin: (np.float32(0.0), np.float32(0.0), np.float32(0.0)) 2025-08-20 15:03:50,331 - Shape (ZYX): (141, 1024, 1024) 2025-08-20 15:03:50,332 - Loading data from outputs/2738_z200to340_segmented_IMM_vertices_normals.csv...
Saved cleaned intensity profiles with features: outputs/2738_z200to340_segmented_OMM_int_profiles_cleaned.pkl Saved statistics: outputs/2738_z200to340_segmented_OMM_filtering_stats.txt
2025-08-20 15:03:50,695 - CUDA GPU is available and will be used 2025-08-20 15:03:50,696 - Starting thickness measurement (1to2)... 2025-08-20 15:03:50,696 - Measuring thickness from surface 1 to surface 2... 2025-08-20 15:03:50,697 - Starting GPU thickness measurement with 314937 points... 2025-08-20 15:03:50,697 - Source points: 139731, Target points: 175206 2025-08-20 15:03:50,698 - Max thickness: 8.0 nm (10.20 voxels) 2025-08-20 15:03:50,698 - Max angle: 1.0 degrees 2025-08-20 15:03:50,706 - Launching CUDA kernel with 1231 blocks and 256 threads per block... 2025-08-20 15:03:50,709 - Retrieving results from GPU... 2025-08-20 15:03:51,038 - Processing matches on CPU to ensure one-to-one matching... 2025-08-20 15:03:54,067 - Found 112416 valid thickness measurements 2025-08-20 15:03:54,069 - Mean thickness: 4.84 nm 2025-08-20 15:03:54,071 - Min: 0.78 nm, Max: 8.00 nm 2025-08-20 15:03:54,072 - Processing completed in 3.38 seconds 2025-08-20 15:03:54,082 - Statistics saved to outputs/2738_z200to340_segmented_IMM_thickness_stats.log 2025-08-20 15:03:54,083 - Creating optimized thickness CSV with only valid measurements... 2025-08-20 15:03:54,130 - Saving 112416 valid thickness measurements to outputs/2738_z200to340_segmented_IMM_thickness.csv 2025-08-20 15:03:56,657 - Thickness statistics: 2025-08-20 15:03:56,659 - Mean: 4.845 nm 2025-08-20 15:03:56,660 - Std: 1.003 nm 2025-08-20 15:03:56,661 - Min: 0.784 nm 2025-08-20 15:03:56,663 - Max: 7.995 nm 2025-08-20 15:03:56,663 - Membrane thickness analysis complete! 2025-08-20 15:03:56,664 - Thickness analysis for IMM completed successfully 2025-08-20 15:03:56,665 - Step 3: Processing intensity profiles for IMM 2025-08-20 15:03:56,666 - Validating segmentation and tomogram compatibility... 2025-08-20 15:03:56,666 - Reading segmentation from inputs/2738_z200to340_segmented.mrc... 2025-08-20 15:03:56,668 - Voxel size: 0.7840 nm 2025-08-20 15:03:56,669 - Origin: (np.float32(0.0), np.float32(0.0), np.float32(0.0)) 2025-08-20 15:03:56,669 - Shape (ZYX): (141, 1024, 1024) 2025-08-20 15:03:56,670 - Reading tomogram from inputs/2738_z200to340.mrc... 2025-08-20 15:03:56,672 - Tomogram voxel size: 0.7840 nm 2025-08-20 15:03:56,672 - Tomogram shape (ZYX): (141, 1024, 1024) 2025-08-20 15:03:56,673 - ✓ Segmentation and tomogram are compatible 2025-08-20 15:03:56,674 - Dimensions: (141, 1024, 1024) 2025-08-20 15:03:56,674 - Voxel size: 0.7840 nm 2025-08-20 15:03:56,675 - ============================================================ 2025-08-20 15:03:56,676 - Processing: 2738_z200to340_segmented_IMM_thickness.csv 2025-08-20 15:03:56,676 - Tomogram: 2738_z200to340.mrc 2025-08-20 15:03:56,677 - ============================================================ 2025-08-20 15:03:56,678 - Loading thickness data... 2025-08-20 15:03:56,847 - Loading tomogram... 2025-08-20 15:03:56,850 - Extracting intensity profiles...
Reading tomogram from inputs/2738_z200to340.mrc... Tomogram voxel size: 0.7840 nm Tomogram shape (ZYX): (141, 1024, 1024) Normalizing tomogram using method: zscore Extracting intensity profiles with extension of 8 voxels Including physical coordinates with voxel size: 0.7840 nm Processing 112416 valid coordinate pairs Extracting intensity values for 112416 profiles...
2025-08-20 15:04:01,646 - Extracted 112416 intensity profiles 2025-08-20 15:04:01,647 - Filtering intensity profiles... 2025-08-20 15:04:01,648 - === Statistics of the extracted intensity profiles ===
Successfully extracted 112416 intensity profiles
2025-08-20 15:04:17,001 - Step 1: Found 48030 candidate profiles 2025-08-20 15:04:17,002 - Median separation: 5.09 ± 1.46 voxels 2025-08-20 15:04:17,003 - Median prominence SNR: 3.7x 2025-08-20 15:04:17,003 - === Apply filtering criteria: (1) two minima (with user-defined SNR compared to baseline) separated by max; (2) the min are positioned between the two matched points from which thickness was measured === 2025-08-20 15:04:37,603 - Step 2: 68312/112416 profiles passed the filtering criteria (60.8%) 2025-08-20 15:04:37,786 - Saving results to: outputs/
=== Summary of intensity profiles post filtering === Total profiles analyzed: 112,416 Candidate profiles found after step 1: 48,030 Profiles that passed criteria: 68,312 (60.8%) Profiles excluded: 44,104 (39.2%) === Characteristics of intensity profiles === Median separation of minima: 5.09 ± 1.46 voxels Median SNR of minima: 3.7x (minima are 3.7x more pronounced than the baseline) === Applied filtering criteria (user-defined) === Required SNR of minima: 0.1x (minima must be 0.1x more pronounced than the baseline) Position filtering: Exact (no margin) Central max required: True Extension range from midpoint: (-8, 8) voxels === Reasons for intensity profile exclusion === minima_outside_extended_region: 35,914 (31.9%) insufficient_snr_min1: 7,814 (7.0%) insufficient_snr_min0: 68 (0.1%) no_central_maximum: 304 (0.3%) zero_baseline_noise: 4 (0.0%) === Quality metrics of included profiles === Mean separation of minima: 5.01 ± 1.34 voxels Median separation of minima: 4.87, Range: 1.85-13.54 Mean SNR of minima: 1249123095870.3 ± 124997031408531.8x Median SNR of minima: 3.0x, Q25-Q75: 2.0-4.2x (Minima are on average 1249123095870.3x more pronounced than the baseline) Saved cleaned thickness DataFrame: outputs/2738_z200to340_segmented_IMM_thickness_cleaned.csv Merging Pass 1 characterization feature data with 112416 profiles Available feature data entries: 112416 Profiles with feature data: 112416/112416 Profiles with minima positions: 48030/112416 Profiles with separation distance: 48030/112416 Profiles with SNR data: 48030/112416 Profiles with central max position: 0/112416 Saved original intensity profiles with features: outputs/2738_z200to340_segmented_IMM_int_profiles.pkl Mapping 68312 cleaned profiles to 68312 passing indices Merging Pass 2 filtering feature data with 68312 profiles Available feature data entries: 112416 Profiles with feature data: 68312/68312 Profiles with minima positions: 68312/68312 Profiles with separation distance: 68312/68312 Profiles with SNR data: 68312/68312 Profiles with central max position: 68312/68312
2025-08-20 15:04:49,881 - Intensity profiling for IMM completed successfully 2025-08-20 15:04:49,883 - Profiles extracted: 112416 2025-08-20 15:04:49,884 - Profiles after filtering: 68312 2025-08-20 15:04:49,886 - Filter pass rate: 60.8% 2025-08-20 15:04:49,887 - Full pipeline completed! 2025-08-20 15:04:49,888 - === Intensity Profiling Summary === 2025-08-20 15:04:49,889 - NE: ✓ Completed - 200594 profiles passed filtering 2025-08-20 15:04:49,889 - OMM: ✓ Completed - 36971 profiles passed filtering 2025-08-20 15:04:49,890 - IMM: ✓ Completed - 68312 profiles passed filtering
Saved cleaned intensity profiles with features: outputs/2738_z200to340_segmented_IMM_int_profiles_cleaned.pkl Saved statistics: outputs/2738_z200to340_segmented_IMM_filtering_stats.txt
In [4]:
# Check results
for membrane_name, result in results.items():
print(f"{membrane_name}:")
print(f" Thickness CSV: {result['thickness_csv']}")
if 'intensity_results' in result:
int_result = result['intensity_results']
print(f" Profiles filtered: {int_result['profiles_filtered']}")
print(f" Pass rate: {int_result['pass_rate']:.1%}")
print()
NE: Thickness CSV: outputs/2738_z200to340_segmented_NE_thickness.csv Profiles filtered: 200594 Pass rate: 90.5% OMM: Thickness CSV: outputs/2738_z200to340_segmented_OMM_thickness.csv Profiles filtered: 36971 Pass rate: 52.2% IMM: Thickness CSV: outputs/2738_z200to340_segmented_IMM_thickness.csv Profiles filtered: 68312 Pass rate: 60.8%
2. CPU-Only Processing¶
For systems without GPU or when CPU parallelization is preferred (parallelization affects only the thickness measurement step; surface extraction and intensity profiling run on single CPU regardless):
In [ ]:
# CPU-parallelized processing
results_cpu = memthick.run_full_pipeline(
segmentation_path=segmentation_path,
output_dir=output_dir,
membrane_labels=membrane_labels,
# Processing options
use_gpu=False,
num_cpu_threads=8, # Adjust for your system
batch_size=1000, # Can be adjusted if memory is an issue but I normally go with 2000
# Surface processing
interpolate=True,
interpolation_points=1,
refine_normals=True,
flip_normals=True,
# Thickness measurement
max_thickness=8.0,
max_angle=1.0,
# Intensity profiling (recommended)
extract_intensity_profiles=True,
tomo_path=tomo_path,
intensity_extension_voxels=8,
intensity_extension_range=(-8, 8),
intensity_min_snr=0.1,
intensity_margin_factor=0.3, # 30% tolerance
intensity_central_max_required=True,
intensity_require_both_minima=True,
# Optional outputs
save_vertices_mrc=False,
save_thickness_mrc=False
)
In [ ]:
surface_results = memthick.process_membrane_segmentation(
segmentation_path=segmentation_path,
output_dir=output_dir,
membrane_labels=membrane_labels,
interpolate=True,
interpolation_points=1,
refine_normals=True,
flip_normals=True,
radius_hit=10.0
)
print("Surface extraction completed:")
for membrane, csv_path in surface_results.items():
print(f" {membrane}: {csv_path}")
Stage 2: Thickness Measurement¶
In [ ]:
# Measure thickness from pre-processed surfaces
for membrane_name in ['NE', 'OMM', 'IMM']:
vertices_csv = f'outputs/2738_z200to340_segmented_{membrane_name}_vertices_normals.csv'
thickness_csv, stats_file = memthick.measure_membrane_thickness(
segmentation_path=segmentation_path,
input_csv=vertices_csv,
output_dir=output_dir,
max_thickness=8.0,
max_angle=1.0,
direction='1to2',
use_gpu=True
)
print(f"{membrane_name} thickness: {thickness_csv}")
Stage 3: Intensity Profile Analysis - Building Intuition¶
- The input tomogram can be normalized using several methods (controlled by
intensity_normalize_method):'zscore': Standardize to zero mean and unit variance (recommended)'minmax': Scale to range [0, 1]'percentile': Clip to 1st-99th percentile range'none': Use original intensity values
- Intensity profiles are extracted along lines connecting paired surface points.
- Profiles are extracted with a pre-set extension range from the measurement midpoint. It's important not to extend them too far, otherwise other structures might be detected (e.g., membrane-bound ribosomes), which could confound peak detection. Controlled by
intensity_extension_voxels(recommended: 8-15) andintensity_extension_range(recommended: (-8, 8) to (-15, 15)). - Profiles can be Gaussian-smoothed if needed to reduce noise (controlled by
intensity_smooth_sigma, where 0=no smoothing). - The detected profile features serve as quality criteria for filtering thickness measurements.
Detected Profile Features¶
- Dual Minima: Represent lipid headgroups (hydrophilic regions) - enable with
intensity_require_both_minima=True - Central Maximum: Represents lipid tails (hydrophobic core) - enable with
intensity_central_max_required=True - Baseline: Calculated from profile edges for SNR estimation - controlled by
intensity_edge_fraction(recommended: 0.1-0.3)
Further Quality Criteria¶
- Signal-to-Noise Ratio: Minima prominence relative to baseline noise - controlled by
intensity_min_snr(recommended: 0.1-2.0) - Minima Positioning: Minima positioned between or near the matched measurement points - controlled by
intensity_margin_factor(recommended: 0.0-0.3, allowing minima to extend slightly beyond measurement points for curved membranes)
In [ ]:
# Apply intensity-based filtering to existing measurements
for membrane_name in ['NE', 'OMM', 'IMM']:
thickness_csv = f'outputs/2738_z200to340_segmented_{membrane_name}_thickness.csv'
filter_results = memthick.int_profiles_extract_clean(
thickness_csv=thickness_csv,
tomo_path=tomo_path,
output_dir=output_dir,
# Profile extraction
intensity_extension_voxels=8,
intensity_extension_range=(-8, ),
intensity_normalize_method='zscore',
# Quality filtering
intensity_min_snr=0.1,
intensity_central_max_required=True,
intensity_require_both_minima=True,
intensity_margin_factor=0.2,
# Smoothing of the profiles
intensity_smooth_sigma=0.0, # Personally, would not smooth
# Outputs
save_cleaned_df=True,
save_profiles=True,
save_statistics=True
)
print(f"{membrane_name} intensity filtering completed")
print(f" Profiles passed: {filter_results['statistics']['profiles_passed']}")
print(f" Pass rate: {filter_results['statistics']['pass_rate']:.1%}")
Command Line Interface¶
Navigate to the /cryocat/ folder and locate memthick.py.
Run Full Pipeline¶
# Full pipeline with intensity profiling (recommended)
python memthick.py \
segmentation.mrc \
--mode full \
--output_dir /path/to/output/
--membrane_labels NE:1,OMM:2,IMM:3 \
--interpolate \
--interpolation_points 1 \
--extract_intensity \
--tomo_path tomogram.mrc \
--max_thickness 8.0 \
--max_angle 1.0 \
--intensity_min_snr 0.5 \
--intensity_margin_factor 0.2 \
--intensity_central_max_required \
--intensity_require_both_minima
Run Separate Pipeline Steps¶
Surface Extraction¶
python memthick.py \
segmentation.mrc \
--mode surface \
--membrane_labels NE:1,ER:2 \
--output_dir /path/to/output/
Thickness Measurement¶
# Requires pre-existing vertices_normals.csv from surface extraction
python memthick.py \
segmentation.mrc \
--mode thickness \
--input_csv vertices_normals.csv \
--max_thickness 8.0 \
--use_gpu
Filter Thickness Results¶
# Requires pre-existing thickness.csv and tomogram
python memthick.py \
segmentation.mrc \
--mode intensity \
--thickness_csv thickness.csv \
--tomo_path tomogram.mrc \
--intensity_min_snr 0.5 \
--intensity_margin_factor 0.2 \
--intensity_central_max_required \
--intensity_require_both_minima
CLI Parameter Reference¶
Required Parameters¶
segmentation_path: Input segmentation file (positional argument)--mode: Pipeline mode (full,surface,thickness,intensity)
Common Options¶
--output_dir: Output directory (default: same as input file)--membrane_labels: Comma-separated label mapping (e.g.,NE:1,ER:2)--config_path: YAML configuration file (alternative to --membrane_labels)--max_thickness: Maximum thickness in nm (default: 8.0)--max_angle: Maximum cone angle in degrees (default: 1.0)--use_gpu/--use_cpu: Processing mode (default: GPU if available)
Intensity Profiling Options¶
--extract_intensity: Enable intensity profiling (flag)--tomo_path: Tomogram file path (required for intensity modes)--intensity_min_snr: Minimum signal-to-noise ratio (default: 1.0)--intensity_margin_factor: Position tolerance, 0.0-0.3 (default: 0.1)--intensity_central_max_required: Require central maximum (flag)--intensity_require_both_minima: Require both minima (flag)
Outputs were uploaded here
In [5]:
NE_thickness_csv = 'outputs/2738_z200to340_segmented_NE_thickness_cleaned.csv'
NE_profiles = 'outputs/2738_z200to340_segmented_NE_int_profiles_cleaned.pkl'
NE_profiles_nonfiltered = 'outputs/2738_z200to340_segmented_NE_int_profiles.pkl'
IMM_thickness_csv = 'outputs/2738_z200to340_segmented_IMM_thickness_cleaned.csv'
IMM_profiles = 'outputs/2738_z200to340_segmented_IMM_int_profiles_cleaned.pkl'
IMM_profiles_nonfiltered = 'outputs/2738_z200to340_segmented_IMM_int_profiles.pkl'
OMM_thickness_csv = 'outputs/2738_z200to340_segmented_OMM_thickness_cleaned.csv'
OMM_profiles = 'outputs/2738_z200to340_segmented_OMM_int_profiles_cleaned.pkl'
OMM_profiles_nonfiltered = 'outputs/2738_z200to340_segmented_OMM_int_profiles.pkl'
In [6]:
# Load thickness files and profiles
NE_data = maple.load_membrane_data(thickness_csv=NE_thickness_csv, intensity_profiles_pkl=NE_profiles)
IMM_data = maple.load_membrane_data(thickness_csv=IMM_thickness_csv, intensity_profiles_pkl=IMM_profiles)
OMM_data = maple.load_membrane_data(thickness_csv=OMM_thickness_csv, intensity_profiles_pkl=OMM_profiles)
Loaded 200594 intensity profiles from outputs/2738_z200to340_segmented_NE_int_profiles_cleaned.pkl Loaded 68312 intensity profiles from outputs/2738_z200to340_segmented_IMM_int_profiles_cleaned.pkl Loaded 36971 intensity profiles from outputs/2738_z200to340_segmented_OMM_int_profiles_cleaned.pkl
In [7]:
# Check loaded data
print("Data Summary:")
for name, data in [("NE", NE_data), ("IMM", IMM_data), ("OMM", OMM_data)]:
summary = data.get_summary()
print(f"{name}: {summary}")
Data Summary:
NE: {'n_thickness_measurements': 200594, 'n_intensity_profiles': 200594, 'thickness_range': (np.float64(1.9204001), np.float64(7.995263)), 'thickness_mean': np.float64(5.715790260377678), 'thickness_std': np.float64(0.6439476034645089), 'coordinate_columns': ['x1_voxel', 'y1_voxel', 'z1_voxel'], 'has_profiles': True}
IMM: {'n_thickness_measurements': 68312, 'n_intensity_profiles': 68312, 'thickness_range': (np.float64(1.7530774), np.float64(7.995263)), 'thickness_mean': np.float64(5.205129820914334), 'thickness_std': np.float64(0.8863175161109151), 'coordinate_columns': ['x1_voxel', 'y1_voxel', 'z1_voxel'], 'has_profiles': True}
OMM: {'n_thickness_measurements': 36971, 'n_intensity_profiles': 36971, 'thickness_range': (np.float64(2.2174869), np.float64(7.995263)), 'thickness_mean': np.float64(4.853210620732467), 'thickness_std': np.float64(0.6645813972741067), 'coordinate_columns': ['x1_voxel', 'y1_voxel', 'z1_voxel'], 'has_profiles': True}
2. Extract Thickness Measurements¶
In [8]:
NE_thickness = maple.analyze_membrane_thickness(
data=NE_data,
outlier_removal_method=None # But can be set to IQR method for filtering the results, percentile as well
)
IMM_thickness = maple.analyze_membrane_thickness(IMM_data)
OMM_thickness = maple.analyze_membrane_thickness(OMM_data)
print(NE_thickness.get_summary())
print(IMM_thickness.get_summary())
print(OMM_thickness.get_summary())
Thickness Analysis Summary: ========================== Number of measurements: 200,594 Mean thickness: 5.72 nm Standard deviation: 0.64 nm Median thickness: 5.71 nm Range: 1.92 - 8.00 nm IQR: 0.86 nm Thickness Analysis Summary: ========================== Number of measurements: 68,312 Mean thickness: 5.21 nm Standard deviation: 0.89 nm Median thickness: 5.02 nm Range: 1.75 - 8.00 nm IQR: 1.26 nm Thickness Analysis Summary: ========================== Number of measurements: 36,971 Mean thickness: 4.85 nm Standard deviation: 0.66 nm Median thickness: 4.83 nm Range: 2.22 - 8.00 nm IQR: 0.71 nm
In [9]:
# Access stats
stats = NE_thickness.statistics
print(f"Mean thickness: {stats['mean']:.2f} ± {stats['std']:.2f} nm")
print(f"Median thickness: {stats['median']:.2f} nm")
print(f"Thickness range: {stats['min']:.2f} - {stats['max']:.2f} nm")
Mean thickness: 5.72 ± 0.64 nm Median thickness: 5.71 nm Thickness range: 1.92 - 8.00 nm
3. Extract Intensity Profile Features¶
In [10]:
# This is the function that works under the hood of the intensity profile plotting function. Not strictly required as standalone
NE_profiles_analysed = maple.analyze_intensity_profiles(
data=NE_data,
extension_range_voxels=(-8, 8), # When you plot the data, this is the distance you'll see from the midpoint
bin_profiles_by_thickness=True, # Bin profiles depending on membrane thickness
thickness_bin_size_nm=0.20 # Bin size
)
IMM_profiles_analysed = maple.analyze_intensity_profiles(
data=IMM_data,
extension_range_voxels=(-8, 8),
bin_profiles_by_thickness=True,
thickness_bin_size_nm=0.20
)
OMM_profiles_analysed = maple.analyze_intensity_profiles(
data=OMM_data,
extension_range_voxels=(-8, 8),
bin_profiles_by_thickness=True,
thickness_bin_size_nm=0.20
)
print(NE_profiles_analysed.get_summary())
print(IMM_profiles_analysed.get_summary())
print(OMM_profiles_analysed.get_summary())
Analyzing 200594 intensity profiles... Created 25 thickness bins for profiles Analyzing 68312 intensity profiles... Created 27 thickness bins for profiles Analyzing 36971 intensity profiles... Created 25 thickness bins for profiles Intensity Profile Analysis Summary: ================================== Number of profiles: 200,594 Extension range: (-8, 8) voxels Interpolation points: 201 Profile type: filtered Binned profiles: Yes Number of thickness bins: 25 Intensity Profile Analysis Summary: ================================== Number of profiles: 68,312 Extension range: (-8, 8) voxels Interpolation points: 201 Profile type: filtered Binned profiles: Yes Number of thickness bins: 27 Intensity Profile Analysis Summary: ================================== Number of profiles: 36,971 Extension range: (-8, 8) voxels Interpolation points: 201 Profile type: filtered Binned profiles: Yes Number of thickness bins: 25
4. Perform a leaflet asymmetry analysis¶
Calculate the ratio of the intensity differences between each minimum and the central maximum (asymmetry score)
Asymmetry scores¶
- Symmetric membranes: Scores near 1.0 (0% asymmetry)
- Asymmetric membranes: Scores e.g., >1.2 (>20% asymmetry)
- Check for outliers: Using the bubble plot or 3D plot
In [11]:
NE_asymmetry = maple.analyze_membrane_asymmetry(
data=NE_data,
thickness_bin_size_nm=0.2, # Bin size for analysis
min_profiles_per_bin=30, # Minimum profiles required for the bin to be considered
use_unity_scoring=True, # 1.0 = symmetric, >1.0 = asymmetric
use_median_aggregation=True,
outlier_removal_method=None, # But outliers can be removed based on IQR filtering
)
print(NE_asymmetry.get_asymmetry_summary())
# Access results
results_df = NE_asymmetry.bin_results
mean_asymmetry = results_df['asymmetry_score'].mean()
print(f"Mean asymmetry score: {mean_asymmetry:.2f}")
IMM_asymmetry = maple.analyze_membrane_asymmetry(
data=IMM_data,
thickness_bin_size_nm=0.2,
min_profiles_per_bin=30,
use_unity_scoring=True,
use_median_aggregation=True
)
OMM_asymmetry = maple.analyze_membrane_asymmetry(
data=OMM_data,
thickness_bin_size_nm=0.2,
min_profiles_per_bin=30,
use_unity_scoring=True,
use_median_aggregation=True
)
Starting asymmetry analysis on 200594 profiles...
Created 31 bins from 1.9 to 8.0 nm
Processed bin 3.42 nm: 55 profiles, asymmetry = 1.326
Processed bin 3.62 nm: 44 profiles, asymmetry = 1.073
Processed bin 3.82 nm: 203 profiles, asymmetry = 1.135
Processed bin 4.02 nm: 478 profiles, asymmetry = 1.050
Processed bin 4.22 nm: 2,021 profiles, asymmetry = 1.065
Processed bin 5.22 nm: 26,371 profiles, asymmetry = 1.069
Processed bin 6.22 nm: 15,461 profiles, asymmetry = 1.058
Processed bin 7.22 nm: 1,885 profiles, asymmetry = 1.036
Successfully calculated asymmetry for 24 bins
Rejected 7 bins (insufficient profiles)
=== Asymmetry Analysis Summary ===
Total profiles analyzed: 200,570
Mean asymmetry: 6.6%
Median asymmetry: 5.5%
Highly asymmetric bins (>20%): 1/24 (4.2%)
Intensity-based calculations: 24
Position-based calculations: 0
Asymmetry Analysis Summary:
==========================
Total bins created: 24
Valid bins: 24
Rejected bins: 7
Total profiles analyzed: 200,570
Asymmetry Statistics:
- Median asymmetry: 1.055 (5.5%)
- Mean asymmetry: 1.066 (6.6%)
- Standard deviation: 0.060
- Range: 1.012 - 1.326
- Notably asymmetric bins (>20%): 1
Mean asymmetry score: 1.07
Starting asymmetry analysis on 68312 profiles...
Created 32 bins from 1.8 to 8.0 nm
Processed bin 2.85 nm: 44 profiles, asymmetry = 1.001
Processed bin 3.25 nm: 206 profiles, asymmetry = 1.173
Processed bin 3.45 nm: 308 profiles, asymmetry = 1.062
Processed bin 3.65 nm: 666 profiles, asymmetry = 1.017
Processed bin 3.85 nm: 2,578 profiles, asymmetry = 1.005
Processed bin 4.85 nm: 3,552 profiles, asymmetry = 1.024
Processed bin 5.85 nm: 3,821 profiles, asymmetry = 1.014
Processed bin 6.85 nm: 1,607 profiles, asymmetry = 1.002
Processed bin 7.85 nm: 339 profiles, asymmetry = 1.015
Successfully calculated asymmetry for 26 bins
Rejected 6 bins (insufficient profiles)
=== Asymmetry Analysis Summary ===
Total profiles analyzed: 68,275
Mean asymmetry: 3.4%
Median asymmetry: 2.3%
Highly asymmetric bins (>20%): 0/26 (0.0%)
Intensity-based calculations: 26
Position-based calculations: 0
Starting asymmetry analysis on 36971 profiles...
Created 29 bins from 2.2 to 8.0 nm
Processed bin 2.92 nm: 39 profiles, asymmetry = 1.159
Processed bin 3.32 nm: 344 profiles, asymmetry = 1.099
Processed bin 3.52 nm: 337 profiles, asymmetry = 1.046
Processed bin 3.72 nm: 310 profiles, asymmetry = 1.040
Processed bin 3.92 nm: 2,674 profiles, asymmetry = 1.022
Processed bin 4.92 nm: 2,231 profiles, asymmetry = 1.005
Processed bin 5.92 nm: 687 profiles, asymmetry = 1.029
Processed bin 6.92 nm: 114 profiles, asymmetry = 1.032
Processed bin 7.92 nm: 47 profiles, asymmetry = 1.062
Successfully calculated asymmetry for 25 bins
Rejected 4 bins (insufficient profiles)
=== Asymmetry Analysis Summary ===
Total profiles analyzed: 36,955
Mean asymmetry: 5.1%
Median asymmetry: 2.9%
Highly asymmetric bins (>20%): 0/25 (0.0%)
Intensity-based calculations: 25
Position-based calculations: 0
In [12]:
thickness_plot = maple.plot_thickness_distribution(
data=[NE_thickness, IMM_thickness, OMM_thickness],
membrane_names=['NE', 'IMM', 'OMM'],
histogram_bins=40,
thickness_range=(1, 10),
show_statistics=True,
show_mean_lines=True,
density_normalization=True,
colors=['#1f77b4', '#ff7f0e', '#2ca02c'],
plot_title='Organelle Thickness Distributions'
)
thickness_plot.show()
thickness_plot.write_html("plots/thickness_distribution.html", include_plotlyjs="cdn")
Statistics for NE: Mean: 5.72 nm Std: 0.64 nm Median: 5.71 nm Count: 200,594 Statistics for IMM: Mean: 5.21 nm Std: 0.89 nm Median: 5.02 nm Count: 68,312 Statistics for OMM: Mean: 4.85 nm Std: 0.66 nm Median: 4.83 nm Count: 36,971
In [ ]:
IFrame("plots/thickness_distribution.html", width=800, height=600)
2. Membrane thickness maps¶
In [16]:
spatial_plot = maple.plot_thickness_3d(
data=[NE_thickness, IMM_thickness, OMM_thickness],
membrane_names=['NE', 'IMM', 'OMM'],
color_scale='OrRd',
color_range=(2, 8), # Thickness color range
marker_size=2,
sample_fraction=0.5, # Downsample for quicker visualization
plot_title='Membrane thickness maps'
)
spatial_plot.show()
spatial_plot.write_html("plots/thickness_3d_plot.html", include_plotlyjs="cdn")
In [ ]:
IFrame("plots/thickness_3d_plot.html", width=800, height=600)
In [ ]:
# Color membranes by their mean thickness
spatial_mean_plot = maple.plot_thickness_3d(
[NE_thickness,IMM_thickness, OMM_thickness],
membrane_names=['NE', 'IMM', 'OMM'],
marker_size=2,
plot_title='Mean membrane thickness map',
color_by_mean=True,
color_range = (4.5,6.5),
color_scale = 'BuPu',
sample_fraction=0.5 # Downsample for quicker visualization
)
spatial_mean_plot.show()
spatial_mean_plot.write_html("plots/thickness_mean_plot.html", include_plotlyjs="cdn")
NE mean thickness: 5.72 nm IMM mean thickness: 5.21 nm OMM mean thickness: 4.85 nm
In [ ]:
IFrame("plots/thickness_mean_plot.html", width=800, height=600)
3. Visualize surface assignment (PCA results)¶
In [17]:
surface_plot = maple.plot_surfaces(
data=NE_thickness,
membrane_names=['NE'],
surface1_color='#1f77b4', # Blue for surface 1
surface2_color='#ff7f0e', # Orange for surface 2
marker_size=3,
sample_fraction=0.5, # Downsample for quicker visualization
plot_title='Assignment of surface points'
)
surface_plot.show()
surface_plot.write_html("plots/surface_assignment_NE.html", include_plotlyjs="cdn")
NE: Surface 1 points: 100,297, Surface 2 points: 100,297
In [ ]:
IFrame("plots/surface_assignment_NE.html", width=800, height=600)
4. Summarized intensity profiles (measurement points are projected as vertical lines)¶
In [18]:
profile_plot = maple.plot_intensity_profile_summary(
data=[NE_profiles_analysed, IMM_profiles_analysed, OMM_profiles_analysed], # Load results of analyze_intensity_profiles(), not directly the .pkl files
membrane_names=['NE', 'IMM', 'OMM'],
extension_range_voxels=(-8, 8), # Generally keep the same as the parameter with which the pipeline was run, otherwise you'll get a straight line at intensity = 0
voxel_size_nm=0.784, # Convert x-axis to nm, otherwise will be shown in voxels
show_individual_profiles=False,
show_mean_profile=True,
show_percentile_bands=True,
colors=['#1f77b4', '#ff7f0e', '#2ca02c'],
plot_title='Extracted intensity profiles (summary)'
)
profile_plot.show()
profile_plot.write_html("plots/intensity_profiles_all.html", include_plotlyjs="cdn")
In [ ]:
IFrame("plots/intensity_profiles_all.html", width=800, height=600)
5. Thickness-binned intensity profiles¶
In [19]:
binned_plot = maple.plot_intensity_profile_binned(
data=NE_profiles_analysed, # Load results of analyze_intensity_profiles(), not directly the .pkl files
membrane_names=['NE'],
thickness_bins=4, # Number of bins
binning_method='quantile', # Equal numbers of profiles per bin
extension_range_voxels=(-8, 8),
voxel_size_nm=0.784,
plot_title='NE profiles per thickness quartile'
)
binned_plot.show()
binned_plot.write_html("plots/intensity_profiles_NE_binned.html", include_plotlyjs="cdn")
Using 4 bins for NE: Bin 1: 1.9-5.3 nm: 56,777 profiles Bin 2: 5.3-5.7 nm: 59,631 profiles Bin 3: 5.7-6.1 nm: 59,829 profiles Bin 4: 6.1-8.0 nm: 54,622 profiles
In [ ]:
IFrame("plots/intensity_profiles_NE_binned.html", width=800, height=600)
In [20]:
equal_width_plot = maple.plot_intensity_profile_binned(
data=NE_profiles_analysed, # Load results of analyze_intensity_profiles(), not directly the .pkl files
membrane_names=['NE'],
thickness_bins=4, # Number of bins
binning_method='equal_width', # Divide the thickness range into bins with equal widths, different number of profiles per bin
extension_range_voxels=(-8, 8),
voxel_size_nm=0.784,
plot_title='NE profiles binned by equal width'
)
equal_width_plot.show()
equal_width_plot.write_html("plots/intensity_profiles_NE_binned_equal_width.html", include_plotlyjs="cdn")
Using 4 bins for NE: Bin 1: 1.9-3.4 nm: 41 profiles Bin 2: 3.4-5.0 nm: 19,851 profiles Bin 3: 5.0-6.5 nm: 160,298 profiles Bin 4: 6.5-8.0 nm: 20,404 profiles
In [ ]:
IFrame("plots/intensity_profiles_NE_binned_equal_width.html", width=800, height=600)
6. Visualize filtering results¶
In [21]:
filtering_plot = maple.plot_excluded_included(
thickness_csv=[NE_thickness_csv, IMM_thickness_csv, OMM_thickness_csv],
membrane_names=['NE', 'IMM', 'OMM'],
filtered_profiles_pkl=[NE_profiles, IMM_profiles, OMM_profiles],
original_profiles_pkl=[NE_profiles_nonfiltered, IMM_profiles_nonfiltered, OMM_profiles_nonfiltered],
plot_title='Excluded vs included point pairs'
)
filtering_plot.show()
filtering_plot.write_html("plots/excluded_vs_included_points.html", include_plotlyjs="cdn")
Loaded 221676 original and 200594 filtered profiles for NE Profiles for NE - Passed: 200594, Failed: 21082 Loaded 112416 original and 68312 filtered profiles for IMM Profiles for IMM - Passed: 68312, Failed: 44104 Loaded 70793 original and 36971 filtered profiles for OMM Profiles for OMM - Passed: 36971, Failed: 33822
In [ ]:
IFrame("plots/excluded_vs_included_points.html", width=800, height=600)
7. Leaflet asymmetry plots¶
In [22]:
# Asymmetry histogram
asymmetry_hist = maple.plot_asymmetry_distribution(
data=[NE_asymmetry, IMM_asymmetry, OMM_asymmetry],
membrane_names=['NE', 'IMM', 'OMM'],
use_percentages=True, # Show as percentages
histogram_bins=30,
asymmetry_range=(0, 30), # Percentage range
y_axis_type='density',
plot_title='Asymmetry distributions'
)
asymmetry_hist.show()
asymmetry_hist.write_html("plots/asymmetry_hist.html", include_plotlyjs="cdn")
In [ ]:
IFrame("plots/asymmetry_hist.html", width=800, height=600)
In [23]:
# Asymmetry vs thickness relationship (to check for outliers)
bubble_plot = maple.plot_asymmetry_vs_thickness_bubble(
data=[NE_asymmetry, IMM_asymmetry, OMM_asymmetry],
membrane_names=['NE', 'IMM', 'OMM'],
asymmetry_range=(0, 30),
plot_title='Asymmetry vs thickness',
)
bubble_plot.show()
bubble_plot.write_html("plots/asymmetry_bubble.html", include_plotlyjs="cdn")
In [ ]:
IFrame("plots/asymmetry_bubble.html", width=800, height=600)
In [24]:
# Spatial asymmetry distribution
asymmetry_3d = maple.plot_asymmetry_3d(
[NE_asymmetry, IMM_asymmetry, OMM_asymmetry],
membrane_names=['NE', 'IMM', 'OMM'],
asymmetry_range=(0,30),
plot_title='Asymmetry 3D distribution',
color_scale='BuPu',
sample_fraction=0.5
)
asymmetry_3d.show()
asymmetry_3d.write_html("plots/asymmetry_3d_plot.html", include_plotlyjs="cdn")
In [ ]:
IFrame("plots/asymmetry_3d_plot.html", width=800, height=600)
8. Plot intensity minima-to-minima distances¶
The minima-to-minima distance is a feature of the profiles that is extracted during analysis and saved in the .pkl file as 'separation_distance'
In [25]:
# Analyze separation between intensity minima
min_min_plot = maple.plot_min_to_min_distribution(
data=[IMM_thickness, OMM_thickness],
membrane_names=['IMM', 'OMM'],
histogram_bins=40,
separation_range=(1, 10), # Distance range
voxel_size_nm=0.784, # Convert to nm, separation_distance is saved in voxels in the pipeline
show_statistics=True,
plot_title='Minima-to-minima distances'
)
min_min_plot.show()
min_min_plot.write_html("plots/min_min_distances.html", include_plotlyjs="cdn")
Statistics for IMM: Mean: 3.92 nm Std: 1.05 nm Median: 3.82 nm Range: 1.45 - 9.65 nm Count: 68,304 Statistics for OMM: Mean: 3.76 nm Std: 0.98 nm Median: 3.66 nm Range: 1.45 - 9.53 nm Count: 36,970
In [ ]:
IFrame("plots/plots/min_min_distances.html", width=800, height=600)
Extract Coordinates of Matched Measurement Points¶
In [26]:
# Extract the coordinates of matched points - can be used to check with the original tomogram
# Returns a dataframe
coord_table = maple.create_profile_coordinates_table(
membrane_data=NE_data,
spatial_bounds={ # Optional, if not specified, will extract all measurement point coordinates
'x': (900, 1000), # X range in voxels
'y': (50, 150), # Y range in voxels
'z': (10, 20) # Z range in voxels
},
include_thickness_data=True,
include_profile_metadata=True,
)
print(f"Found {len(coord_table)} profiles in region")
coord_table.head()
Found 1558 profiles in region
Out[26]:
| profile_index | p1_x | p1_y | p1_z | p2_x | p2_y | p2_z | thickness_nm | passes_filter | has_features | profile_length | minima_between_points | separation_distance | prominence_snr | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 9727 | 976 | 54 | 10 | 979 | 60 | 8 | 5.488000 | True | True | 24 | True | 4.173913 | 3.553441 |
| 1 | 9730 | 971 | 56 | 10 | 975 | 62 | 7 | 6.123236 | True | True | 25 | True | 4.666667 | 8.587264 |
| 2 | 9733 | 972 | 56 | 10 | 974 | 62 | 7 | 5.488000 | True | True | 24 | True | 5.565217 | 4.650548 |
| 3 | 9736 | 973 | 56 | 10 | 974 | 62 | 8 | 5.020050 | True | True | 24 | True | 4.869565 | 3.060674 |
| 4 | 9738 | 966 | 59 | 10 | 969 | 65 | 7 | 5.761200 | True | True | 25 | True | 5.333333 | 7.102457 |
Convert Thickness Data Into Motive Lists¶
Can be useful for ChimeraX visualization
In [27]:
maple.save_thickness_motls(NE_thickness_csv, output_directory='outputs/')
maple.save_thickness_motls(IMM_thickness_csv, output_directory='outputs/')
maple.save_thickness_motls(OMM_thickness_csv, output_directory='outputs/')
Saved motls with 200594 points Saved motls with 68312 points Saved motls with 36971 points